微生物共发生网络中节点度的幂律分布和在R中拟合（更正）

微生物共发生网络中节点度的幂律分布和在R中拟合

非常抱歉先前的文档出了计算错误，修改后重新发送下，还望大家见谅。（感谢大佬指出问题，先前错在置换检验那一步中，以此篇为准）

在提到微生物共发生网络分析的文献中，通常会在附录中出现类似这样的图，绘制了网络节点度的频数分布，并拟合了幂律分布，告知大家网络图的节点度具有这种分布特征。

微生物共发生网络中节点度的幂律分布

一般而言，微生物共发生网络的度符合幂律分布，即大部分物种具有少量的连接数，极个别的物种具有非常多的连接数，是典型的无标度网络（Steele et al. 2011），说明微生物群落构建方式是非随机的（Barberan et al. 2012）。

在构建了微生物共发生网络之后，需要首先检查网络是否符合这种分布。如果服从，说明构建的网络无问题，符合一般规律，可以进行后续的统计分析；如果拒绝，则此时需要小心，是计算过程出了问题，选择的方法不合适，还是数据来自于特殊的情况（上述是一般规律，但不否认非常规的存在）。

R语言拟合节点度的幂律分布

接下来就展示在R中统计微生物共发生网络中，节点幂律分布的方法。

示例数据和R代码链接（提取码：nys3）：

https://pan.baidu.com/s/1TuA9iNXe8VLzZS4hfu6aTg 

准备文件

首先需要提供一个记录网络中节点度频数分布的文件。

#读取度频数分布文件
dat <- read.delim('degree.txt')

* 节点度频数分布计算

如果还没有计算，则应首先计算网络中各节点的度，并获得度频数分布统计，节点度的概念及计算可见“节点和边特征”。

网络文件以邻接矩阵作为输入，文件基本格式及igraph包的入门操作可见前文网络基础概述。

#通过 igraph 包计算获得节点的度频数分布
library(igraph)

#输入数据示例，邻接矩阵
#这是一个微生物互作网络，数值“1”表示微生物 OTU 之间存在互作，“0”表示无互作
adjacency_unweight <- read.delim('adjacency_unweight.txt', row.names = 1, sep = '\t', check.names = FALSE)
head(adjacency_unweight)[1:6] #邻接矩阵类型的网络文件

#邻接矩阵 -> igraph 的邻接列表，获得非含权的无向网络
igraph <- graph_from_adjacency_matrix(as.matrix(adjacency_unweight), mode = 'undirected', weighted = NULL, diag = FALSE)
igraph #igraph 的邻接列表

#计算节点度
V(igraph)$degree <- degree(igraph)
head(V(igraph)$degree)

#度分布统计
degree_dist <- table(V(igraph)$degree)
degree_num <- as.numeric(names(degree_dist))
degree_count <- as.numeric(degree_dist)

dat <- data.frame(degree = degree_num, count = degree_count)
head(dat)

#查看度分布
#可观察到微生物相关网络通常服从幂律分布
par(mfrow = c(1, 3))
hist(V(igraph)$degree, xlab = 'Degree', ylab = 'Frequency',
main = 'Degree distribution')
plot(degree_count, degree_num, xlab = 'Degree', ylab = 'Count',
main = 'Degree distribution')
plot(degree_count, degree_num, log = 'xy', xlab = 'Log-degree',
ylab = 'Log-count', main = 'Log-log degree distribution')

通过节点度的频数分布拟合幂律分布

在R中，可通过非线性回归函数nls()拟合幂律分布。

#拟合，a 和 b 的初始值手动指定，数据不同需要多加尝试
mod <- nls(count ~ a*degree^b, data = dat, start = list(a = 2, b = 1.5))
summary(mod)

#提取关键值
a <- round(coef(mod)[1], 3)
b <- round(coef(mod)[2], 3)
a; b

可知幂律方程式为“y = a*x^b = 224.754x^-1.325”。

不过nls()函数的拟合结果不带有R²，因此R²还需额外计算。可以使用构建好的方程，通过预测变量（这里为degree）的值拟合响应变量（这里为count）的值，再根据预测值和观测值的差异，获得R²。

#使用构建好的方程，通过预测变量拟合响应变量，再根据和观测值的差异，获得 R2
#SSre：the sum of the squares of the distances of the points from the fit
#SStot：the sum of the squares of the distances of the points from a horizontal line through the mean of all Y values
fit <- fitted(mod)
SSre <- sum((dat$count-fit)^2)
SStot <- sum((dat$count-mean(dat$count))^2)
R2 <- round(1 - SSre/SStot, 3)
R2

显示R²是很高的。

那么，该方程式“显著”吗？同样地，由于nls()函数的拟合结果不带有P值，因此也需要额外计算。

这里就可以通过“万能”的置换检验，计算P值。

#p 值可以根据置换检验的原理获得
#将 count 的值分别随机置换 N 次（例如 999 次），通过随机置换数据后数据获取 R2（R2'）
#比较随机置换后值的 R2' 大于观测值的 R2 的频率，即为 p 值
p_num <- 1
dat_rand <- dat
for (i in 1:999) {
dat_rand$count <- sample(dat_rand$count)
SSre_rand <- sum((dat_rand$count-fit)^2)
SStot_rand <- sum((dat_rand$count-mean(dat_rand$count))^2)
R2_rand <- 1 - SSre_rand/SStot_rand
if (R2_rand > R2) p_num <- p_num + 1
}
p_value <- p_num / (999+1)
p_value

P值小于0.01，是很显著的。

最后可知，该微生物共发生网络的节点度分布符合幂律分布，方程式“y = a*x^b = 224.754x^-1.325”，R²=0.983，P<0.001。

ggplot2绘制节点度分布图

最后使用ggplot2，将度分布更形象地绘制出，并添加拟合线，以及将计算结果标注在图中。

#作图展示
library(ggplot2)

p <- ggplot(dat, aes(x = degree, y = count)) +
geom_point(color = 'blue') +
theme(panel.grid.major = element_line(color = 'gray'), panel.background = element_rect(color = 'black', fill = 'transparent')) +
stat_smooth(method = 'nls', formula = y ~ a*x^b, method.args = list(start = list(a = 2, b = 1.5)), se = FALSE) +
labs(x = 'Degree', y = 'Count')

#添加公式拟合的注释
label <- data.frame(
formula = sprintf('italic(Y) == %.3f*italic(X)^%.3f', a, b),
R2 = sprintf('italic(R^2) == %.3f', R2),
p_value = sprintf('italic(P) < %.3f', p_value)
)

p + geom_text(x = 13, y = 210, aes(label = formula), data = label, parse = TRUE, hjust = 0) +
geom_text(x = 13, y = 190, aes(label = R2), data = label, parse = TRUE, hjust = 0) +
geom_text(x = 13, y = 170, aes(label = p_value), data = label, parse = TRUE, hjust = 0)

多网络的分面图作图

最后再展示一个组合样式。

####一个多网络的分面图
#数据
dat2 <- read.csv('degree2.csv', stringsAsFactors = FALSE)

#构建计算函数
lm_labels <- function(dat) {
#拟合
mod <- nls(count ~ a*degree^b, data = dat, start = list(a = 2, b = 1.5))
a <- round(coef(mod)[1], 3)
b <- round(coef(mod)[2], 3)

#计算R2
fit <- fitted(mod)
SSre <- sum((dat$count-fit)^2)
SStot <- sum((dat$count-mean(dat$count))^2)
R2 <- round(1 - SSre/SStot, 3)

#计算 p 值（置换检验原理）
p_num <- 1
dat_rand <- dat
for (i in 1:999) {
dat_rand$count <- sample(dat_rand$count)
SSre_rand <- sum((dat_rand$count-fit)^2)
SStot_rand <- sum((dat_rand$count-mean(dat_rand$count))^2)
R2_rand <- 1 - SSre_rand/SStot_rand
if (R2_rand > R2) p_num <- p_num + 1
}
p_value <- p_num / (999+1)

#方程式值的列表
data.frame(formula = sprintf('italic(Y) == %.3f*italic(X)^%.3f', a, b),
R2 = sprintf('italic(R^2) == %.3f', R2),
p_value = sprintf('italic(P) < %.3f', p_value))
}

#计算及作图
library(plyr)
library(ggplot2)

label <- ddply(dat2, 'Type', lm_labels)

p <- ggplot(dat2, aes(x = degree, y = count, color=Type)) +
geom_point() +
facet_wrap(~Type, nrow = 1) +
stat_smooth(method = 'nls', formula = y ~ a*x^b, method.args = list(start = list(a = 2, b = 1.5)), se = FALSE, show.legend = FALSE) +
labs(x = 'Degree', y = 'Count')

p + geom_text(x = 40, y = 180, aes(label = formula), data = label, parse = TRUE, hjust = 0, color = 'black', show.legend = FALSE) +
geom_text(x = 40, y = 160, aes(label = R2), data = label, parse = TRUE, hjust = 0, color = 'black', show.legend = FALSE) +
geom_text(x = 40, y = 140, aes(label = p_value), data = label, parse = TRUE, hjust = 0, color = 'black', show.legend = FALSE)

参考资料

Barberán, Albert, Bates S T, Casamayor E O , et al. Using network analysis to explore co-occurrence patterns in soil microbial communities. The ISME Journal, 2012, 6(2):343-351.

Steele J A, Countway P D, Xia L, et al. Marine bacterial, archaeal and protistan association networks reveal ecological linkages. The ISME Journal, 2011, 5(9):1414-1425. 

友情链接

相关性和网络分析基础

网络模块内连通度（Zi）和模块间连通度（Pi）及在R中计算 

降维分析

非约束排序（描述性的探索性分析）：

非约束排序中被动添加解释变量：被动添加解释变量

约束排序（将解释变量通过回归方程拟合响应变量的统计模型）：

RDA、CCA的变差分解（VAP） 

对称分析（这类方法意在描述两个或多个矩阵之间的相关性）：

多元因子分析（MFA）

监督降维（带监督的降维方法，常用于分类）：

判别分析（DA）

聚类和分类

层次聚类（无监督，描述性的探索性分析）：

层次聚合：层次聚合分类    层次聚类结果的比较和评估

层次分划：双向指示种分析（TWINSPAN） 

非层次聚类（无监督，描述性的探索性分析）：

划分聚类：k均值划分（k-means）    围绕中心点划分（PAM）

模糊聚类：模糊c均值聚类（FCM）

避免不存在的类

潜变量分类（无监督，潜变量也可视为某种意义上的“降维”）：

潜类别分析（LCA）    潜剖面分析（LPA）约束聚类（无监督，将解释变量通过回归方程“约束”响应变量的模型）：多元回归树 

监督分类（通过已知先验分组构建分类器模型）：

决策树    随机森林分类